Mathematical modeling of intratumoral immunotherapy yields strategies to improve the treatment outcomes

Intratumoral injection of immunotherapy aims to maximize its activity within the tumor. However, cytokines are cleared via tumor vessels and escape from the tumor periphery into the host-tissue, reducing efficacy and causing toxicity. Thus, understanding the determinants of the tumor and immune response to intratumoral immunotherapy should lead to better treatment outcomes. In this study, we developed a mechanistic mathematical model to determine the efficacy of intratumorally-injected conjugated-cytokines, accounting for properties of the tumor microenvironment and the conjugated-cytokines. The model explicitly incorporates i) the tumor vascular density and permeability and the tumor hydraulic conductivity, ii) conjugated-cytokines size and binding affinity as well as their clearance via the blood vessels and the surrounding tissue, and iii) immune cells—cancer cells interactions. Model simulations show how the properties of the tumor and of the conjugated-cytokines determine treatment outcomes and how selection of proper parameters can optimize therapy. A high tumor tissue hydraulic permeability allows for the uniform distribution of the cytokines into the tumor, whereas uniform tumor perfusion is required for sufficient access and activation of immune cells. The permeability of the tumor vessels affects the blood clearance of the cytokines and optimal values depend on the size of the conjugates. A size >5 nm in radius was found to be optimal, whereas the binding of conjugates should be high enough to prevent clearance from the tumor into the surrounding tissue. In conclusion, development of strategies to improve vessel perfusion and tissue hydraulic conductivity by reprogramming the microenvironment along with optimal design of conjugated-cytokines can enhance intratumoral immunotherapy.


Interstitial pressure-fluid velocity
Normal and tumor tissues have properties similar to those of a porous medium.Brinkman's equation describes the flow in porous medium in cases where the velocity gradients are non-negligible and reads: Where µ is the dynamic viscosity, ρ is the density, p i is the interstitial pressure, v f is the fluid velocity and k th is the hydraulic conductivity of the interstitial space [24].For incompressible fluid flow the conservation of mass reads: Where the first term of Q m describes the fluid flux entering from the blood vessels and the second term the flux exiting through the lymphatic system.L p is the blood vessels' hydraulic conductivity, and p v is the vascular pressure.L pl , S vl and p vl are the corresponding parameters for the lymphatic vessels [23].
Intratumorally injected conjugated-cytokines The free conjugated-cytokines that travel in the tumor interstitial space, Ic f , can be transferred due to convection and diffusion, where D Ic f is the diffusion coefficient of the drug in the interstitial space and v f is the fluid velocity.Moreover, the free conjugated-cytokines are transferred across the tumor blood vessel and lymphatic vessel wall (Q Ic ).The remaining terms describe the binding and unbinding of the conjugated-cytokines: c e is the concentration of collagen, k on , k of f are the binding and unbinding rate constants, respectively and Φ is the volume fraction [14].

Conjugated-cytokines in the blood compartment
The injected conjugated-cytokines concentration in the blood is described as, Where V −Q Ic dV dom is the total rate of mass of the agonist that can be transported from the tumor and host tissue to the blood, V blood is the volume of blood and δ clear the rate of clearance from blood.

Transport from tumor and host tissue to blood
Transport across the vessel and lymphatic vessel wall based on Starling's approximation: where P er is the vascular permeability of the conjugated-cytokines, σ f the reflection coefficient and δ Ic f the rate constant that describes the rate in which the conjugated-cytokines exit through the lymphatic vessels.The parameters L p , P er and σ f are expressed as a function of the vessel wall pores and the size of the conjugated-cytokines [7,17]: where γ is the fraction of the vessel wall surface area occupied by pores, r 0 the pore radius, µ the viscosity and L vw the thickness of the vessel wall.H and W describe the steric and hydrodynamic interactions of the conjugated-cytokines with the pores of the vessel wall that hinder diffusive and convective transport respectively and D 0 is the diffusion coefficient of a particle in free solution given by the Stokes-Einstein equation.By ignoring electrostatic interactions H and W become [7]: where F is the partition coefficient expressed as: where λ is the ratio of the conjugated-cytokines size to the vessel wall pore size and K t and K s are expressed as [7]:

Molecular radius of agonist
The radius of the injected agonist is expressed as [22]: where M w is the molecular weight of agonist.

Boundary conditions
Needle outlet boundary conditions Where V in is the inflow velocity, Ic f In is the initial concentration of injected cytokines.
Needle periphery boundary conditions External surface boundary conditions Modeling of immune response and tumor growth

Implementation of model equations
The model equations are solved using COMSOL Multiphysics (COMSOL, Inc.Burlington, MA, USA).At the tumor center, the initial cancer cell density was assumed to have its peak value and then it decreases alongside the straight line due to a step function.The domain above a threshold cancer cell density is assumed as the tumor region and below that, the host tissue.
The parameters change from their abnormal-tumor value to their normal value as a function of the cancer cell density.The values change at the threshold cancer cell density using a step function.The step functions are used to certify continuity of the model variables.
Equations were solved at a 1D geometry with spherical symmetry.At the left side of the 1D interval is the tumor center.As we move away from the tumor center we move towards the host tissue.Furthermore, previously used equations 6-15 are also applied on the immune response and tumor growth model.

Kinematics of tumor growth
The growth stretch ratio is calculated as [2]: Where λ g is the growth stretch ratio R T the rate of change of the concentration of cancer cells and T 0 the initial concentration of cancer cells.The tumor growth is implemented using deformation of the spatial frame mesh relative to the material frame mesh with a prescribed mesh displacement of Where x the spatial frame coordinates and X the material frame coordinates.The solid velocity is calculated as: Where t is the time.

Interstitial pressure-fluid velocity
Normal and tumor tissues have properties similar to those of a porous medium.
According to Darcy's law and the mesh movement due to solid velocity, the interstitial fluid velocity is given by: where k th is the hydraulic conductivity of the interstitial space [24].The mass balance gives [27], [4], The first term of the right-hand side of the equation describes the fluid flux entering from the blood vessels and the second term the flux exiting through the lymphatic system.L p is the blood vessels' hydraulic conductivity, and p v is the vascular pressure.L pl , S vl and p vl are the corresponding parameters for the lymphatic vessels [23].
Intratumorally injected conjugated-cytokines The free conjugated-cytokines that travel in the tumor interstitial space, Ic f , can be transferred due to convection and diffusion, where D Ic f is the diffusion coefficient of the conjugated-cytokines in the interstitial space and v f is the fluid velocity.Moreover, the free conjugated-cytokines are transferred across the tumor blood vessel and lymphatic vessel wall (Q Ic ).The remaining terms describe the binding and unbinding of the conjugated-cytokines: c e is the concentration of collagen, k on , k of f are the binding and unbinding rate constants, respectively and Φ is the volume fraction [14].

Pro-inflammatory cytokines from immune cells
The pro inflammatory cytokines produced by the immune cells can be transported by convection and diffusion: Where the right-hand side terms describe the production of pro inflammatory cytokines by innate immune cells, effector CD8+ and CD4+ Tcells and antigen presenting cells.
The last term describes the degradation of cytokines.

Total pro-inflammatory cytokines
The total pro-inflammatory cytokines are a combination of the injected cytokines, the injected cytokines bound to collagen and the pro-inflammatory cytokines produced by the immune cells.
The cytokines related to the injected agonist are solved in mol m 3 and the cytokines by the immune cells in kg m 3 thus the molecular wight M w converts the units to enable the summation.

Trafficking of immune cells
Normalization of the tumor micro-environment increases trafficking of immune cells when combined with immunotherapy.When the vascular density doubles from 50cm −1 to 100cm −1 the source of immune cells increases by 1.367.By assuming that the function is 1 for the host tissue value of (70cm −1 ).The trafficking functions reads pro-inflammatory cytokines produced by the immune cells [18].
where S v represents the functional vascular density measured in cm −1 .

5/19
Immature antigen presenting cells The immature antigen presenting cells are expressed as: Where the first right-hand term describes the source of immature antigen presenting cells, the second term the degradation and the last term the reduction due to activation.The activation to antigen presenting cells depends on the pro-inflammatory cytokines and the interaction of immature antigen presenting cells with the tumor cells and antigen.

Antigen presenting cells
The antigen presenting cells are expressed as: Where the first right-hand term describes the increase of antigen presenting cells due to the activation from immature antigen presenting cells and the last term is a degradation term.

Effector CD4+ T cells
The effector CD4+ T cells are expressed as: Where the first term describes the source of effector CD4+ T cells which is assumed to be analog to the concentration of antigen presenting cells responsible for the activation of CD4+ T cells in the lymph nodes.The last term is the degradation of CD4+ T cells.

Effector CD8+ T cells
The effector CD8+ T cells are expressed as: Where the first term describes the source of effector CD8+ T cells which is assumed to be analog to the concentration of antigen presenting cells responsible for the activation of CD8+ T cells in the lymph nodes.The last term is the degradation of CD8+ T cells.

Innate cells
The innate immune cells that induce cytolysis are expressed as: Where the first term describes production of innate cells which depends on the concentration of pro-inflammatory cytokines.

Cancer cells
The concentration of cancer cells is expressed as: Where the first right-hand side term describes the proliferation of cancer cells due to oxygen and the last term describes the killing of cancer cells by innate cells, immature antigen presenting cells and effector CD8+ T cells.

Antigen
The concentration of cancer cells is expressed as: Where the first right-hand side term describes the production of antigen induced by the cytolytic effect of innate cells and effector CD8+ T cells.The last term describes the antigen uptake by the immature antigen presenting cells.

Oxygen transport
The rate of change of oxygen concentration in the tissue was modeled with a convection diffusion equation that includes a source and a sink term [20], [11].The source term is due to oxygen supply from the blood vessels and the sink term describes oxygen consumption by cancer cells: where S v is the vascular density, D ox the oxygen diffusion coefficient, A ox and k ox are oxygen uptake parameters, c iox is the oxygen concentration in the vessels, v f is the fluid velocity and P erox is the vascular permeability of oxygen defined as the oxygen diffusion coefficient divided by the length of the vessel wall.

Immune cell death rate
According to experimental data [3], a 40 times decrease in oxygen concentration (from 20% to 0.5%) doubled the apoptotic rate of immune cells.Thus the degradation rates are expressed as: Hydraulic conductivity (host) Oxygen concentration in the vessels Oxygen diffusion coefficient Oxygen uptake parameter A Oxygen uptake parameter k

Table A .
Mathematical model characteristics compared to other models

Table B .
Table of model variables

Table C .
Table of model parameters